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ABSTRACT 


The  Alfven  problem  has  plagued  explicit  numerical  calculations 
of  magnetohydrodynamic  fluid  behavior  in  regions  where  the  fluid 
density  is  relatively  small  but  can  be  solved  by  the  inclusion  of 
terms  representing  the  displacement  current  in  Maxwell's  Equations. 
This  modification  limits  the  important  physical  velocities  to  c  in  a 
realistic  way  and  co  permits  much  longer  computational  timesteps  than 
had  previously  been  possible.  Thus  the  small  c  technique  will  greatly 
extend  the  usefulness  of  present  day  computers  for  solving  these  MHD 
problems . 


PROBLEM  STATUS 
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PHYSICALLY  MOTIVATED  SOLUTION  OF  THE  ALFV3&  PROBLEM 


J.  P,  Boris 

Naval  Research  Laboratory 
Washington,  D.C.,  20390 

The  Alfveh  problem  has  plagued  explicit  numerical  calculations  of  magneto¬ 
hydrodynamic  fluid  behavior  in  regions  where  the  fluid  density  is  relatively 
small  but  can  be  solved  by  the  inclusion  of  terms  representing  the  displace¬ 
ment.  current  in  Maxwell's  Equations,  This  modification  limits  the  important 
physical,  velocities  to  c  in  a  realistic  way  and  so  permits  much  longer  compu¬ 
tational  time-steps  than  had  previously  been  possible.  Thus  the  small  c 
technique  will  greatly  extend  the  usefulness  of  present  day  computers  for 
solving  these  MHD  problems. 


For  years,  the  Alfven  problem  has  seriously  limited  explicit,  multidimen¬ 
sional  computations  of  the  magnetohydrodynamic  (MHD)  behavior  of  conducting 
fluids  where  the  mass  density  p  or  the  magnetic  field  strength  B  varies  greatly 
from  one  point  in  the  calculations!  region  to  another.  Problems  arise  because 

i 

the  Alfven  velocity, 

—A  E  2  /  'nr^'  (1) 

can  vary  widely  across  a  highly  inhomogeneous  magneto-fluid  system  requiring  a 

prohibitively  small  compute ional  tiaestep,  5t,  to  provide  numerical  stability 

1  2 

according  to  the  usual  stability  criterion,  * 

6t  <  6x/max  (VA(x,t)),  (2) 

where  the  maximum  is  taken  over  the  whole  system.  This  stability  condition 
s bates  in  general  that  numerical  information  about  changes  in  the  fluid  must 
propagate  faster  than  the  largest  of  the  characteristic  fluid  velocities.  In 
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the  Alfven  problem  imder  consideration  h?re,  is  the  largest  velocity  of  interest 
principally  because  p  becanes  very  small  in  certain  regions,  say  at  the  edges 
of  a  magnetically  confined  plasma  or  at  the  top  of  a  gravitationally  stratified 
atmosphere.  A  small  time step  is  not  intrinsically  disadvantageous,  of  course, 
since  phenomena  happening  on  the  dynamic  timescales  of  the  Alfven  waves  in  the 
low  density  regions  may  be  of  interest.  Usually,  however,  the  physics  of  pri¬ 
mary  interest  occurs  in  the  high  density  regions  where  most  of  the  fluid  is  lo¬ 
cated.  If  the  central  density  is  four  orders  of  magnitude  greater  than  the  den¬ 
sity  of  the  fluid  at  the  vail,  as  can  occur  in  seme  plasma-pinch  calculation 
for  instance,  one  may  be  faced  with  a  timest^p  which  is  100  times  smaller  than 
desjred  simply  to  provide  numerical  stability  in  a  region  of  negligible  physical 
interest  anyway. 

The  stability  condition  (2)  applies  to  explicit  calculations  where  numerical 
"information"  can  propagate  only  to  the  r.earest-neighbor  cells  during  a  single 
timestep.  This  defines  a  numerical  velocity 

V  =  6x/6t. 
num 

Thus  a  numerical  instability  in  an  explicit  computational  model  is  much  like  a 
shock.  The  inability  of  the  computational  fluid  to  propagate  energy  away  frem 
an  accumulation  point  due  to  a  slow  numerical  speed  means  that,  a  perturbation  in 
&he  ccanputed  flow  can  accumulate  unendingly  until  the  calculation  must  be  stopped. 
By  increasing  the  numerical  velocity  until  it  is  greater  than  all  the  characteris¬ 
tic  velocities  of  the  fluid,  this  numerical  "steepening"  of  perturbations  can  be 
removed.  This  usually  involves  shortening  6t  however. 

This  Alfven  problem  is  confronted  in  all  explicit  numerical  calculations  of 


strongly  inhomogeneous  magnetofluids  when  a  partial  differential  equation  describ¬ 
ing  the  evolution  of  the  magnetic  field  must  be  solved  simuiirneously  with 
equations  describing  the  fluid  motion.  Various  methods  of  circumventing  the 
problem  have  been  attempted  previously  but  none  has  been  fully  satisfactory. 
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Using  an  implicit  rather  than  explicit  finite  -  difference  formulation  of  the 

5  4 

relevant  equations  can  have  the  effect  of  raising  to  infinity,  5  but  usual¬ 
ly  becar.es  prohibitively  complicated  in  more  than  one  dimension.  The  program¬ 
ming  almost  invariably  becomes  difficult  and  iterations  may  be  necessary,  hence, 
computational,  effort  increases  greatly  per  tiroestep(even  though  the  time3teps 
can  be  taken  over  a  somewhat  longer  physical  time). 

In  explicit  calculations  various  artificial  corrections  have  been  used  but 
the3e  lack  a  physical  basis  and  so  the  properties  of  the  fluid  are  often  grossly 
altered.  For  instance,  artificial  limits  on  the  density  set  an  upper  limit  to 
the  Alfven  velocity  but  seriously  hinder  fluid  motions  in  the  artificially  dense 
regions,  can  cause  anomalous  instabilities  of  a  convective  type,  and  may  damage 
conservation  of  density.  Another  less  than  ideal  fix  involves  applying  an  impli¬ 
cit  numerical  diffusion  term  to  the  explicit  difference  scheme.  The  propagation 
speed  of  wave  energy  due  to  the  numerics  can  still  be  too  small  for  stability,  but 
the  implicit  diffusion  will  prevent  instabilii  r  by  diffusing  away  perturbation 
accumulations  before  they  can  become  dangerously  large.  Unfortunately,  this  has 
the  effect  of  changing  wavelike  phenomena  to  evanescent  or  diffusion  waves  wher¬ 
ever  the  numerical  speed  becomes  too  small. 

In  nature ,  the  Alfven  problem  is  solved  by  limiting  all  important  velocities 
to  c,  the  velocity  of  light.  For  the  Alfven  and  magnetosonic  waves  of  interest 
here  this  is  done  by  the  displacement  current  in  Maxwell ' s  equations  which  is 
usually  neglected  in  MUD  models.  By  including  this  term,  an  improved  MHD  model 
can  be  f  und  which  does  not  suffer  the  Alfven  problem.  Consider  the  following 
equations : 


ft  = 


dV  = 


dt 


7  P  + 


J.  x  £. 


d_ 

dt 


(f/py)  =  0 


(5) 
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1  5JB.  _ 
c  dt 


V  X  E 


J  _  I  _  _  1  3E 

7  -T^l*  I  -lS77t  . 


(3)  (con't) 


Using  the  definition  of  the  current  J,  the  continuity  aquation,  and  Ohm's  law 
in  the  simple  form  E  +  V  x  B  =  0.  The  momentum  anri  magnetic  equations  can  be 


rewritten  as 


=  V_  x  (V  x  B) 


where  J>  I  (pT  +  -g-)S  +  PV.V  -  -j& 


is  the  usual  pressure-stress  tensor.  The  magnetic  equation  is  unchanged.  The 
last  term  in  the  momentum  equation  (^a)  arises  from  the  displacement  current.  This 
term  may  be  rearranged  using  the  magnetic  equation  to  give  a  modified  momentum 


equation, 


where 

*  SS?  [<*  *  £)A  x  S)3-  l5LPr’  %,] 

is  the  displacement-current  correction  to  the  stress  tensor. 

The  physical  importance  of  Eq.  (5)  arises  from  the  inclusion  of  the  elec¬ 
tromagnetic  effects  in  conservative  form.  The  time  rate  of  change  of  a  generalized 
momentum  is  equal  to  the  divergence  of  a  generalized  pressure-stress  tensor.  While 
the  correction  2P ^  is  essentially  relativistic ,  being  of  order  (v/c)2  relative  to 
ether  stress  tensor  terms,  the  Poynting  vector  contribution  to  the  generalized 
momentum  is  basically  nonrelativistic,  depending  only  on  large  fields  or  low 
density  to  be  important.  In  fact,  if  the  electromagnetic  field  is  thought  of  as 
having  mass,  one  car  write  a  generalized  mass-density  tensor  £*  frem  Eq.  (5), 


x  _EL)  x 
If-rc2 


=  -V  •  -  V_  •  IPj. 


b 


(6) 


R2  BiB  i 

d  *  =  o  6  +  u  6  -  9  • 

ij  ij  id  4rrca 

The  generalized  momentum  is  -.hen  _p*  •  V  • 

Hie  computational  importance  cf  Ea.  (5)  can  be  seen  best  by  solving 
for  the  linear  dispersion  relations  of  the  Alfven  and  magnetosonic  model 
which  are  found  from  the  following  linear  equation5 

-  +  w2  (v  XV.)  X  V,  +  (S2+V  2)(k*V)k 


Here  S2  s  is  the  square  of  the  sound  speed  and  (Eq.  (l))  is  the  vector 
Alfven  velocity.  Assuming  that  k  is  perpendicular  to  both  V  and  £  gives  the 
dispersion  relation  for  transverse  Alfven  waves, 

(f) » =  K  ^  (8) 

where  Y  =  (1  +  /  c2)  1  and  where  0  is  the  angle  between  It  and  B.  As  p  approaches 

zero,  Y  approaches  c2/V^  and 

(f)  lA  **  c2cos265  (8a) 

a  finite  rather  than  infinite  phase  velocity. 

The  four  magnetosonic  (magnetcacouetic )  roots  are  found  under  the  assump¬ 
tion  V_  lies  in  the  plane  established  by  k  and  B.  After  seme  algebra 

KS  54  (s2cos20  +  YS2sin20  +  "fJp/2 

X^/'( S2cca2  +  YS2sin20  +  YV|)2  -4ycos20S2v|. 

This  reduces  to  the  usual  result' ’  when  Y  =  1  and,  when  y  becomes  small 
(VA/c  »  1), 


(h*4)  ^*4)1  "  (4*V)k  -  (k*V}V^j=  0. 
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FMS 

2 

SMS 


-  c2, 


c2Sacos29 

ca+S2eos2fi 


(10) 


where  F  and  S  refer  to  "Fast"  and  "Slow11  magnetosonic  modes  respectively.  Notice 
that  the  slow  mode  is  primarily  ccmpressional  parallel  to  _B  and  does  not  propagate 
across  the  magnetic  field.  As  for  the  Alfven  modes,  all  velocities  are  less  than 
or  equal  to  c. 

Various  levels  of  numerical  approximation  are  possible.  The  Alfven  problem 
will  be  solved,  of  course,  if  Eqs,  (4),  (5)  and  (6)  are  used  as  they  stand,  but 
the  relativistic  correction  is  not  strictly  necessary  unless  (v/c)2  ~  1  in 
which  case  the  definition  of  fluid  momentum  must  be  slightly  modified.  Thus  JP^ 
can  usually  be  dropped  to  greatly  simplify  computation.  Another  simplification 
can  be  made  which  further  reduces  calculation.  The  off  diagonal  term  B^B_.  in 
2*  can  be  set  to  zero.  This  reduces  calculation  greatly  but  gives  linearized 
dispersion  relations  which  are  not  strictly  correct.  The  slow  magnetosonic 
modes  are  primarily  acoustic  (compressional)  when  ^  S.  Since  the  use  of  a 
scalar  P*  slows  down  ail  modes  equally  in  low  density  regions,  these  compressional 
modes  are  slowed  but  really  should  not  be.  This  effect* although  spurious,  is 
not  likely  to  be  very  important  in 

(1)  Low  density  regions  where  an  artificially  small  value  of  c  may  be 
required  to  allow  sufficiently  long  tiraesteps  and  where  the  physics  is  often 
of  little  interest  anyway; 

(2)  High  density  regions  wher®  ~  S,  where  c  is  much  larger  than  either, 
and  thus  the  difference  of  p*  from  p  is  negligible  anyway. 

The  idea  of  usir.g  an  artificially  small  value  of  c  is  well  established, 
being  rather  like  the  use  of  mass  ratio  100  for  raanyboay  plasma  simulations  in¬ 
volving  the  inf i”idual  trajectories  of  many  ions  and  electrons.  In  the  present 
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case  all  of  the  important  physical  effects  are  preserved  and  can  be  scaled 
analytically  to  the  correct  value  of  c  in  those  problems  where  an  artificial 
value  is  used  for  c. 

Figure  1  shows  the  numerical  Alfven-wave  phase  velocity  (x's)  computed 

/  7 

using  the  scalar  p*  correction  for  a  simple  plane-wave  Alfven  excitation.  The 

phase  velocities  relative  to  c  are  plotted  versus  V^/c;  the  solid  curve  is  the 

theoretical  dispersion  relation.  The  asymptote  is  also  shown.  In 

the  configuration  (B  =  constant,  P  =  constant,  V  perturbation  sinusoidal  in  X) 

x  y 

the  dispersion  relations  are  as  given  earlier  for  the  transverse  Alfven  and  fast 
magnetos cnic  modes.  The  slow  magnetosonic  modes  will  of  course  be  slowed  arti- 
fically  since  the  off  diagonal  g*  terms  have  been  neglected. 

Figure  2  shows  the  perturbation  amplitude  as  a  function  of  time.  The  slow 
decay  of  the  wave  shows  the  effect  of  a  small  stabilizing  viscosity.  The  fre¬ 
quencies  and  hence  phase  velocities  shown  in  Fig.  1  were  taken  frctn  a  series  of 
graphs  similar  to  Fig.  2.  The  systematic  difference  between  the  numerical  points 
and  theoretical  curve  in  Fig.  1  is  of  order  3$  and  is  easily  explained  by  the 
finite  difference  approximations  employed  for  the  temporal  and  spatial  sinusoids. 

A  physical  solution  to  the  Alfven  problem  in  explicit,  multidimensional 
computer  calculations  of  low  density  MHD  fluids  has  been  proposed  and  tested. 

The  inclusion  of  terms  arising  from  the  physical  displacement  current  in  Maxwell’s 
equations  limits  physically  important  velocities  to  c,  the  velocity  of  light,  and 
therefore  much  larger  computational  timesteps  can  be  used  than  was  heretofore 

possible.  Factors  of  10  have  been  realized  in  calculations  where  p  /  P  .  ~  104. 

max  min 


Acknowledgement 

The  author  would  like  to  thank  Dr,  J.  M.  Dawson  of  Princeton  University  and 


7 


Marshal  Greenblatt  of  the  Naval  Research  Laboratory  for  many  useful  suggestions 
and  interesting  discussions. 

References 

"hc.V,  Roberts  and  D.E.  Potter,  Methods  in  Computational  Physics,  Vol.  9  (Eds. 

Alder,  Fernbach,  and  Rotenberg,  Academic  Press,  New  York,  1970). 

2 

R.  Courant,  K.O.  Friedrichs,  and  H.  Lewy,  Math.  Ann.  100,  (1923). 

^R.D.  Richtmyer  and  K.W.  Morton,  Difference  Methods  for  Initial  Value  Problems 
(interscience  Publishers,  New  York,  1967). 

^A.  Ralston,  A  First  Course  in  Numerical  Analysis,  (McGraw  -  Hill  Book  Company, 

New  York,  1965). 

'  J.D.  Jackson,  Classical  Electrodynamics  (John  Wiley  &  Pons,  Inc.,  New  York, 

1962). 

^A.G.  Kulikovsky  and  C.A.  Lyubimov,  Magi.etchydrod/namics  (Addison  -  Wesley 
Publishing  Company,  Inc.,  Palo  Alto,  1965). 

7 

K.V.  Roberts  and  J.P.  Boris,  Trinity:  Programs  for  5b  Magnetohydrodynamics , 
in  Proceedings  of  the  Conference  on  Computational  Physics ,  Culham  Laboratory, 

July  28-30,  1969.  Tiie  calculations  reported  here  were  performed  by  the 
MR  HYDE  program,  developed  by  the  author  at  Princeton  and  Culham. 

Computer  Credit 

The  computations  reported  on  here  were  performed  on  the  IBM  360/91  computers 
at  Princeton  University  and  the  Johns  Hopkins  Applied  Physics  Laboratory.  The 
former  computer  facilities  were  supported  in  part  by  the  National  Science  Foundation 
Grant  NSF-GP  579. 


8 


ALFVEN  WAVE  DISPERSION  RELATION 
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The  systematic  error  between  theoiy  (solid  line)  and  computation  (X’s)is 
accounted  for  by  the  finite-difference  representation  of  trigonometric 
functions. 


0  8  !6  24  32  40  48  56  64  72 

t - ► 


Fig.  2  -  Computational  Alfven-wave  amplitude  as  a  function  of  time.  Non¬ 
zero  viscosity  and  thermal  conductivity  account  for  the  slow  decay  of  the 
oscillation. 
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